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Nature of energy transport around a critical point is studied in the Creutz cellular automaton. 
. Fourier heat law is confirmed to hold in this model by a direct measurement of heat flow under a 

temperature gradient. The thermal conductivity is carefully investigated near the phase transition 
by the use of the Kubo formula. As the result, the thermal conductivity is found to take a finite 
value at the critical point contrary to some previous works. Equal-time correlation of the heat flow 
is also analyzed by a mean-field type approximation to investigate the temperature dependence of 
thermal conductivity. A variant of the Creutz cellular automaton called the Q2R is also investigated 
and similar results are obtained. 



CO 



a 

C 

o 
o 



oo 

On 



X 



05.60.+w,44.10.+i,05.50.+q 



I. INTRODUCTION 



Creutz devised a deterministic dynamics for the two-dimensional Ising model with a momentum term Q]. This 
dynamics is a kind of cellular automaton (CA), where states are upadated in a deterministic way with energy con- 
servation and we call it the Creutz cellular automaton (CCA). In the CCA random numbers are not necessary for its 
time evolution, which provides an advantage in numerical simulations. Thus, the CCA and its variants have been used 
t~H , to investigate equilibrium properties of magnetic systems || instead of conventional Monte Carlo method, especially 
t-H ■ the critical phenomena at the critical point. 

Besides this advantage, the CCA provides a dynamical model for time evolution with energy conservation. Thus 



the CCA can be used to study transport phenomena where flows of physical quantities take important roles. In fact, 
numerical results for heat conduction in the CCA were reported in M and the thermal conductivity was found to 



be proportional to T~ 2 in high temperature, where T denotes the temperature. Harris and Grant showed that this 
temperature dependence is explained by the Kubo formula J3|. They presented an asymptotic expression for the 
thermal conductivity in the high and low temperature limits by evaluating the first term in the Kubo formula. 

Then, it is natural to ask a possible connection between the thermal conductivity and the phase transition. Because 
the specific heat diverges at the critical point of the Ising model, the thermal conductivity may also show some 
peculiarity at the point. Actually, in some materials, abnormal behavior of the thermal conductivity has been observed 
||. Clearly, the CCA is suitable to look into the thermal conductivity near the critical point. In the above mentioned 
paper, Harris and Grant made a comment that the thermal conductivity must vanish at the critical temperature Tq 
without any evidences. It is a purpose of our paper to clarify what really happens to the energy transport at the 
critical point in the CCA. 

Thus, in this paper we investigate the temperature dependence of the thermal conductivity in the thermodynamic 
limit. We obtain the thermal conductivity by two methods. One is a direct measurement of heat flow under a 
temperature gradient. The validity of the Fourier heat law is established in a wide range of temperature values and 
the coefficient of thermal conductivity is estimated. The other is the use of the Kubo formula. Explicit derivation of 
the formula is given and the coefficient of thermal conductivity is calculated from equilibrium autocorrelations of the 
energy flow. We check that both the methods yield the same result and find a finite conductivity at To, which does 
not agree with the previous belief. 

We also develop a mean-field approximation for the equal-time correlation of the energy flow, which improves the 
estimate by Harris and Grant 0]. Since it is the first term of the Kubo formula, the result of this treatment not 
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only explains the temperature dependence at the high and low temperature limits, but also gives a qualitatively good 
estimate for the overall temperature dependence. 

The conditions under which the Fourier heat law is satisfied have been studied in the literature mainly by 

using Hamiltonian systems. The dynamical rules of CA are so simple and local that fast simulations are possible. 
Thus, one of the present authors (ST) applied CA to this problem and found some rules that clearly satisfy the Fourier 
heat law However, most of the studies have so far been confined in one-dimensional models, which might cause 
pathological effects due to a single path of the flow. Here we study a two dimensional system with CCA where we 
are free from the above anxiety. We find that the Fourier heat law holds at all temperatures in the CCA. We also 
find that the Q2R, which is a variant of the CCR without momentum terms, satisfies the Fourier heat law in two 
dimensions, although energy transport in Q2R is ballistic in one dimension. 

This paper is organized as follows. In Sec. 2, our model and method are explained and an expression for the local 
energy flux is derived via the equation of continuity. In Sec. 3, we demonstrate that the Fourier heat law holds in a 
wide range of temperatures by a direct simulation. In Sec. 4, the thermal conductivity is calculated by the use of the 
Kubo formula and its temperature dependence is carefully investigated especially around Tq. A mean field analysis 
is done for the correlation of the energy flux in Sec. 5. Numerical results for the Q2R are exhibited in Sec. 6. We give 
summary and discussion in Sec. 7. 



II. MODEL 



The CCA is defined as follows. Let us consider the square lattice. A couple of variables (<7jj, d - ^-) are assigned 
at a site Here Uij G {+1,-1} denotes a spin and cr^j S {0,1,2,3} is called a momentum. Then the total 

Hamiltonian is given by 

id id 

The first term represents the ferromagnetic Ising interaction between the nearest-neighbor spins and the second term 
represents a kind of kinetic energy. We divide the lattice into two sublattices like the checkerboard. Site is called 
even or odd according as the sum i + j is even or odd. One unit of time evolution consists of two processes each of 
which simultaneously updates variables on a sublattice. Namely, when the variables are updated from the states at 
time t, first the even sites are updated at time t + 1/2 and next the odd sites are updated at time t+1. The updating 
rule is the following. Spin flip is accepted when the momentum at the site can compensate the energy change of the 
flip. That is, if the following relation is satisfied, < cr-jj — \<Jij Ylnn ann — ^' wnere nn denotes the nearest neighbor- 
sites of the spin (jjj changes its sign and the momentum is changed to conserve the total energy. 



Now we derive expressions for a local energy and an energy flux. From the total Hamiltonian (2.1), we can define 
the local energy on the site at time t as 

E if j — —0i,j{<Ji+\,j + + (2.2) 

Note that the total energy is equal to the sum of the local energies over the lattice. First we consider the case where 
the site is even. If the spin at site (i,j) is flipped at time t + 1/2, we have 

CT *j -~ a idi 

~ t + 1/2 ~t I f f t , Jt I _t I _t \ 

a i,J = a id ~ ^ ' + + *>J — 1 + a 'i,j + l)' 

and the difference between the local energies at times t + 1/2 and t is given by 

E% 1/2 -Etd=-KM-id+<i-J- 
If the spin crjj is not flipped, the local energy does not change. Thus the energy change is generally expressed as 

p t+l/2 pi „ t , t t \ r/ t+V 2 i _t \ 

\j ~ \j = - 2cr ' i jl cr -i-lj + a i,j-l) d i a id + a hj> 

t /t + 1/2 t \ , * / t+1/2 t \ fr, o\ 

where S(x) is Kronecker's delta 
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S(x) 



1 if i = 
otherwise 



(2.4) 



and we have used the equality 8{x + y) = (1 — xy)/2 that holds for {+1, —1}- The energy difference between 

t + 1/2 and t + 1 is calculated in the same manner, and we obtain 



E 



t+i 



pt+l/2 t+1/2 t+1/2 s , t+1 , t+l/2s . t+1/2 t+1/2 r/ t+1 



t+1/2 



l) 



t+1/2/ t+1/2 



t+1 > 
r i+l,jJ 



- (T 



t+1/2/ t+1/2 



h3 



\"i,j+l "i,j + l> 



(2.5) 



Combining Eqs. (2.2) and (2.5), we obtain the following expression for the energy difference between t and t + 1. 



Elt+l 



,t+i 



where we have used the fact that a 



I _t ^t+1 \ I ^.t + 1/^.t ,_t+l \ 

t+1/2 _ „t r_ nAA /, ' \ _ , t+1/2 



(2.6) 



erf ^ for odd (i, j) and <r id 
energy is conserved, Eq. (2.6) must represent the equation of continuity, 
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er*^ 1 for even Because the total 



E i,3 



Tpt jt 



J, 



t,3,x 



J, 



i,3 + l,V 



h3,y> 



(2.7) 



wher e J- ■ a (a — x or y) denotes the a component of the energy flux at site at time t. Comparing Eqs. ( |2.6| ) 
and (2.7), we find that the components of the energy flux are given as 





-^-lj^j 
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_ t+i i t+i 
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i,j+l,v 


_ t+i, t+i 


- ^ij + l)- 



The same argument can also be applied to the case where site is odd. As the result we arrive at the following 
expressions for the energy flux. If site (i,j) is even, 



jt _ t („t + i _ _t ^ 

"i,j,x — a i-l,j\ a i.j a i.j) 

jt _ _t t-t+1 _ t \ 
J i.i,y — a i,i-l\ a i.i "U! 



(2.8) 
(2.9) 



and if site is odd, 



T t _ t+i i -t+i _ t \ 

J i,3,x — a i—l,j\ a i,j a i,jn 

jt _ t+1 (t+1 _ t \ 

J h3,y ~ a i,j-l\ a i,j a i,3>- 



(2.10) 
(2.11) 



III. THERMAL CONDUCTION UNDER A GRADIENT OF THE TEMPERATURE 

In this section, we report numerical results on energy transport in the CCA obtained by a direct simulation. We 
took the systems of size L x L where L varies from 10 to 300. The periodic boundary condition was imposed on 
the y direction. At the ends in the x direction, two heat reservoirs, one at temperature X"l and the other at Tr, 
were attached as shown in Fig. 1. Each heat reservoir consisted of spins on two vertical lines, where the spins on 
a sublattice were simultaneously updated by the use of the Monte-Carlo method with the heat-bath algorithm. We 
have numerically confirmed that if the two heat reservoirs have an identical temperature the system relaxes to the 
equilibrium state at that temperature. This relaxation to the equilibrium was also observed in the case where only 
one reservoir was attached to the system. 

Energy transport occurs when the left and right reservoirs have different temperatures. It is found that the relaxation 
time to a stationary state is very long at low temperatures below Tq = 2/ log(l + y/2) ~ 2.270, while it is rather short 
at high temperatures. 

The following two cases are examined with particular care. One is the case where both the temperatures Tl and 
T R are higher than T c , T L = 5.0 and T R = 5.5. This is called case A. The other is the case of T L = 2.1 and T R = 2.2, 
where both the temperatures are lower than To We call this case B. In each case, within 10 7 time steps the system 
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of any size (L < 300) reached a stationary state where a uniform flux in the x direction is realized. After the system 
reached the stationary state, we continued the simulation by 10 7 more steps for which we took time averages of 
physical quantities. 

First we consider the distribution of a local kinetic energy, Pij (a) . Because the system is translation invariant in 
the vertical direction, we computed the average of Pij over the vertical line and found that it is given by a canonical 
distribution 



1 L 

3=1 



where Pi is a fitting parameter which is regarded as the local inverse temperature at line i. Figures 2(a) and (b) 



show the distributions for case A and case B, respectively. They clearly demonstrate the property (3.1). Thus local 
equilibrium is realized and the local temperatures are well defined. 

Let Ti denote the temperature at horizontal position i, namely Tj = /?r . We plotted as a function of x = i/L 
for various Ls in Fig. 3(a) and (b), which correspond to cases A and B, respectively. Clearly the scaling limit 

T(x) = lim T\ Lx ], (3.2) 

L — >oc 

where [Lx] means the integer part of Lx, exists and is smooth in both the cases A and B. 
Next we observed the total energy flux per row in the stationary state 



L 



1,3- 



where Jtot,x is the total energy flux in the x direction and the bars mean the time average in the stationary state. If 
the Fourier heat law is realized, this quantity must converge to a nonzero constant in the limit L — > oo with Tl and 
Tr fixed, because then this quantity is written as 



Ji 



tot, a; 



n{T)dT (3.4) 

T L 

with use of the thermal conductivity k(T). We utilized this property to judge whether the Fourier heat law is satisfied 
or not. 

In Fig. 4, the L dependence of Jtot.x/L is shown for various temperature values. The figure shows that the size 
dependence disappears in the large systems. Thus we conclude that the Fourier heat law is realized in a wide range 
of temperatures including the critical point. 

Moreover, Jtot,x/L has a finite value and changes smoothly around the critical temperature. This observation 
suggests that the thermal conductivity has no strong singularity at Tq. However, we can treat not the thermal 
conductivity itself but the integration of it in the present method and a possible singularity, if any, is hardly observed. 
Thus in the next section we investigate the thermal conductivity in the bulk at a given temperature with use of the 
Kubo formula. 



IV. THERMAL CONDUCTIVITY COMPUTED VIA THE KUBO FORMULA 



According to the Kubo formula, the thermal conductivity is equal to the summation of the equilibrium autocorre- 
lation functions of the energy flux as 

1 °° / 1 \ 

<T) = E^tV^tot J 1 - 2 5 *.o . (4-1) 
t=o ^ ' 

where J\ ot ^ x = Y^,ij ^i,j,x 1S the total energy flux in the x direction at time t, (• • •) means the equilibrium ensemble 
average at temperature T, and N is the total number of sites. This formula is proved for the CCA in Appendices A 
and B. 

We numerically computed the autocorrelation functions (J t ° ot x J* ot x ) for t < 150 in the CCA under the periodic 
boundary conditions in the x and y directions. Initial conditions were randomly generated by a Monte-Carlo method 
with temperature T. We denote the partial Kubo sum up to time t by «;*, namely 
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Kt = £> ( J ™,A X ) (i - b v , Q ) . (4.2) 

t'=0 ^ ' 

Figure 5 shows numerically computed «;* in the system of size 200 x 200 at various temperatures. It is observed 
that the summation converges by t = 30 for every temperature. At temperatures above Tq, the sum monotonically 
increases and tends to a constant exponentially fast. At low temperatures the monotonicity is lost and significant 
fluctuations appear. However, k* still reaches a convergence by t — 10. 

Figure 6 shows the thermal conductivity thus obtained and that computed via the direct measurement of energy 
flux as explained in the previous section. Both the results agree with each other very well. From the figure we know 
that the thermal conductivity has a peak at T ~ 2.70, which is slightly above the critical temperature Tq. Above 
the peak value, the thermal conductivity gradually decreases and tends to zero in the high temperature limit. Below 
the peak value, the conductivity shows a remarkable change around Tc and reaches nearly zero at T = 2.0. Detailed 
measurements were done near the critical temperature Tq and the results are shown in Fig. 7. This figure shows that 
k(T) appears continuous and smooth at the critical point, though the magnitude of the change is large. Because little 
size dependence is seen when L > 100, we can conclude that at least no divergence or no vanishment of k(T) occurs 
at To Of course we cannot deny the possibility of singularity or discontinuity in a higher derivative. 



V. MEAN-FIELD ANALYSIS OF THERMAL CONDUCTIVITY 



In this section, we estimate the equal time correlation function of the heat flow using a mean-field approximation 
and discuss its temperature dependence. This quantity is the first term of the Kubo formula Eq. (4.1), namely k°(T), 
and thus we can obtain some information on the temperature dependence of the thermal conductivity. 

As derived in Appendix B, k°(T) is expressed in terms of an average of the total flow Jtot,x in the local equilibrium 

as 

K ° (T) = ^v< J tot,*-W> = lT J% H0 N(n-T R y (5 - 1} 

where (• • -)i c denotes the average with respect to the local equilibrium product measure ( |A6| ) with the right reservoir 
temperature Tr, and the left reservoir temperature Tl- 

First we consider the quantity (Jj,,-, x )i e at an even site Denoting the local equilibrium measure by p\ e , we 

have 

{Ji,j,x)\c ^ ^ JiJ.xPlc- (5-2) 
{*.*} 



Substituting Eq. (2.8) into Jij, x , we can express this as 

(Ji,j,x)le = X! Irr < l -..'' 7 ',; - °'i-l,j a 'i,j)Ple 
{*,*} 

E* 
,j <Ti,j Pie, (5.3) 

{<?,*} 

where a[ a denotes the updated spin value at the even site and J^* means the summation over the configurations 
in which the spin aij can flip. Whether the spin flip occurs or not depends on the spin and the momentum variable 
at and the sum of the spin values on the nearest neighbor sites, 

h = <7i_ij + (r i+ ij + (Ji.j-i + (5.4) 

Specifically, the spin flip is possible in the following configurations: 

(±4, ±1,(3, 2)) 
(±2, ±1,(3, 2,1)) 

{h,<T iJ ,o id )={ (±0, ±1,(3, 2, 1,0)) (5.5) 
(±4,T-1,(1,0)) 
(±2,T-1,(2,1,0)) . 
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Because the summation (5.3) must be taken over the configurations for the whole system, it is difficult to be carried 
out exactly. Thus we consider the following mean field approximation. In this approximation the spin variables at 
the next nearest neighbor sites are replaced by their average values. Those average values should depend only on the 
horizontal position and not on the vertical position, since the local equilibrium measure is translation invariant in 
the y direction. Thus the average concerning the local equilibrium measure is replaced by the average concerning the 
following measure 



PjOjV, (<t)) 

E*,*p(w ( cr ))' 

where 



(5.6) 



with 



P(ct,<t, (cr)) = exp(-4/3iCTij) exp[a<7i_i j + ba hj -i + ca hj+i + da i+1 ,j 

+ "'i.iCA'-i^-iJ + + ftoij-i + PiVi+tj)] (5.7) 



a = f3i- 2 (<j)i-2 + A-i(cr)i_i + /3i_i(<7)j_i, 

b = j3i-i{a)i-x + Pi(<j)i + 0i(a) i+ i, 
c = Pi-i(a)i-i + Pi(<j)i + 0i(a) i+ i, 
d = (3 i+ i(a) l+1 +(3 i+ i(a) i+ i + p i+1 (a) i+2 . 



The summation in the denominator of ( |5.6| ) is taken over possible values of and the nearest neighbor spins. 

Here Pi is the inverse temperature at horizontal position i and P[_i takes the same value as ft-i. We introduced P\_i 
for later convenience. (a)i denotes the local equilibrium value of the spin variable at horizontal position i. 
Under the above approximation, (Jij. x )\c is represented as 

2 dZ* 

{Ji,j,x)i e - —g -Qpi ~ > ( 5 - 8 ) 

where Z and Z* are defined by 

a, a 

Z = ^P( C r,a,(a)). 

<7,<7 

With straightforward calculation Z* is obtained as, 

Z* = 2e /3 -i-' 3 * (e~ 4ft + e~ m ) cosh(a + b + c + d) 

+ 2(e- im + e" 6ft + e~ 2ft ) jeft-' 3 -! C osh(a - b - c - d) + e /3 - 1 " ft cosh(a -b + c + d) 

+e^- 1_ft cosh(a + b - c + d) + e^-i~ ft cosh(a + b + c- d)\ 
+ 2(e~ 12ft +e- 8 ^ +e- 4 ^ + 1) 

x {e '- 1 ^ cosh(a + b - c - d) + e /3 *- 1 ~' 3i cosh(a - b - c + d) + e /3 '- 1 " /3i cosh(a -b + c-d) 

+ e ~^~ Pi cosh(-a - b + c + d) + e+^~ Pi cosh(-a + b - c + d) + e _/3 *~ ft cosh(-a + b + c - d)\ 

+ 2e -^-i+A (e~ 4ft + e- m ) cosh(a + b + c + d) 

+ 2(e~ 10ft + e~ 6ft + e- 213 ') je^-i-ft cosh(-a + b + c + d) + c~^-^ cosh(a -b + c + d) 
+e -Pi-i+Pi cosh(a + b-c + d) + e-P'i-x+Pi C osh(a + b + c - d)} . (5.9) 



In the first order of AT(:= T, - T;_i) 



dZ* 



/3' 1=/ 3 i +/3 2 AT can De simplified as 



G 



dZ* 



— rp2 

/3;_ 1 =/3 i + / 3^AT 



4AT 

{( e -4ft +e -8ft )cosh(12/3i((j)i) 



+ 4(e- 10ft + e- 6ft + e" 2A ) cosh(6A(^)i) + 3(1 + e~ 8ft )(l + e- 4ft )}- (5.10) 

Z is also calculated as 

Z = 2 4 (cosh 4 (3ft (ct), + ft) + cosh 4 (3ft (a), - #)) (1 + e - 4ft )(l + e~ 8ft )- (5.11) 
Thus we arrive at the approximate formula for k (T), 

k°(T)= lim (Jtot f /iV)le 

cosh(12/3(cr)) cosh(2/3(cr)) + 2(1 + 2 cosh 4/3) cosh(6/3(cr)) + 6 cosh 2/3 cosh 4/? 



4T 2 [cosh 4 (3/3(cr) +/3) + cosh 4 (3/3(cr) - /?)] cosh 2/3 cosh 4/3 



(5-12) 



where in the limit of AT — * the system becomes uniform and we identify (Jtot.x/N)\ c with (Jij_ x )ic and put 
(a) := (a)i. In the high temperature limit, using (a) = 0, we have 



~ ( 5 - 13 ) 



while in the low temperature limit, using (a) = 1, 



k°(T) ~ ± e - 8 /*\ (5.14) 



These asymptotic forms are th e same as obtained by Harris and Grant 



However, the formula (5.12) gives more information about overall temperature dependence. Although the present 
approximation is not good near the critical point, within this approximation we find that K (T) is continuous but 
shows a cusp at the mean- field critical temperature Tm — 3.5 because (a) on (T M -T)2. In Fig. 8 we compare the 
mean-field results with the numerical ones obtained in the previous section. In the high temperature region both the 
results agree with each other, while discrepancies appear at low temperatures. This is partly due to the difference 
between the mean-filed critical temperature and the true critical temperature Tq — 2.27. In addition the simulation 
results have no cusp and actually change smoothly. 

In Fig. 9 we show k(T) and 3.5 x k (T) both of which are numerically obtained from the equilibrium autocorrelation 
functions of the energy flux. This figure shows that k°{T) is nearly proportional to k(T) in high temperatures. This 
implies that the autocorrelation functions of the energy flux are similar in this temperature region, which is also 
perceived by comparing the two curves for T = 3.0 and T = 3.5 in Fig. 5. 



VI. HEAT CONDUCTION IN THE Q2R 

As a simplified variant of the CCA, the Q2R was devised and some equilibrium and dynamical features were 
investigated |ll],[l^,|l0| . There are no momentum variables in the Q2R, where a spin flips only when the sum of the 
nearest-neighbor spins is zero. Despite the simplicity of the model, it is known that the critical behavior for the 
magnetization can be simulated by this model. 

We have performed direct simulations of the Q2R in contact with two heat reservoirs at different temperatures in 
almost the same manner as in Sec. 3. Heat reservoirs were realized by the same algorithm as shown in Fig. 1. The 
temperatures of the reservoirs were set as Tl = 6.0 and Tr = 10.0. Here we took quasi one-dimensional systems of 
s ize L x 1 wit h various Ls. Simulation time for each size is 5 x 10 7 time steps. The expressions for the energy flux 
and (|2~10| ) can be used without changes because they do not contain momentum variables. For the same reason, 



in the Q2R we cannot determine local temperature from the distribution of local kinetic energy as was done in the 
CCA. Thus we plotted local energies in the stationary state for various system sizes in Fig. 10. As in the Creutz 



model, the Q2R also shows a smooth energy profile in the scaling limit (3.2). System-size dependence of the total 
energy flux is shown in Fig. 11. The total energy flux converges to a nonzero finite value in the limit L — » oo and it 
demonstrates that the Q2R has a normal thermal conductivity at least when the temperatures are sufficiently high. 
This means that the normal thermal conductivity in the CCA is not caused by the presence of the momentum terms. 
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The thermal conductivity was carefully calculated with use of the Kubo formula in a system of size 100 x 100 at 
temperatures around Tq. The result is shown in Fig. 12, which exhibits similar behavior to the CCA. The thermal 
conductivity shows a remarkable change near Tq but seems continuous and smooth. This result disagrees with Costa 
and Herrmann [^o| , where they reported that energy flux vanished at the critical point and no transport occurred 
below the critical point. This discrepancy may be attributed to the differences in system sizes and heat reservoirs in 
their and our systems. In [flCf] , the distance between the reservoirs is 20, which may be too small to obtain bulk thermal 
conductivity. Their heat reservoir is deterministic and keeps energy a constant in the boundary layer representing 
the reservoir. Thus the motion of the total system must eventually turn into periodic. Because energy flow rarely 
occurs in low temperature, such simple dynamics may not be able to generate it, whereas our reservoirs are stochastic 
and rare events can happen. Another possible interpretation is that they misunderstand the great change of thermal 
conductivity around the critical point as vanishing. 

In addition, Costa and Herrmann reported two different types of transport processes. One is normal diffusion and 
the other is a systematic transport called "highway." The latter causes a ballistic transport. However, we did not 
find such ballistic transport in our simulations. This is also attributed to the differences in heat reservoirs and system 
sizes. The highway is characteristic of their deterministic reservoirs and moreover the fraction of highways decreases 
to zero as the system size increases. 

At the end of this section, we mention the one-dimensional Q2R dynamics. If i is even, the spin value of site i at 
time t + 1 is expressed in terms of spin variables at time t as 



= a! +1/2 = at 



2ct S i 



cr,_i a, a. 



i+l! 



'i+U 



(6.1) 



In the same manner, if i is odd 



t+ l {+1/2 {+1/2 x ( {+1/2 . t+l/2\ 



Defining local energy of site i at time t as 



_{ _{ JtJt t 



we obtain the following relation for the local energy using Eqs.(6.1) and (|6.2[). Namely if i is even, 



crf+l _ rpt 



(6.2) 
(6.3) 
(6.4) 



and if i is odd, 



(6.5) 



Therefore, the energy transport in one-dimensional Q2R is ballistic and the Fourier heat law is not satisfied. Thus 
we have found that the dimensionality has an important role for the Fourier heat law in the Q2R. 



VII. SUMMARY AND DISCUSSION 



In this paper we have studied the thermal conduction in the CCA with two methods. One is the direct measurement 
of the heat flux under a temperature gradient. The other is the use of the Kubo formula. The former revealed that 
the assumption of local equilibrium is satisfied and that Fourier heat law is realized in a wide range of temperatures. 
The thermal conductivity was carefully calculated near the critical point by the latter method and the results show 
no singularity for k(T) at To 

How a thermal conductivity behaves at Tq is a highly nontrivial problem. Harris and Grant || and Costa and 
Herrmann flpfl both argued that the thermal conductivity vanishes at the critical point. On the other hand, the 
autocorrelation of the total energy flux might show a slow decay due to the critical slowing down. Then the thermal 
conductivity might be divergent at To Our present result shows either is not the case. 

The present result does not mean that there is no singularity in energy transport at Tq ■ The Fourier heat law means 
that the macroscopic motion of energy density obeys the diffusion equation with diffusion constant D(T) = k(T) /C(T), 
where C(T) is the specific heat. The present result shows that k(Tq) is finite while C(T) diverges to infinity at To 
Thus the diffusion constant D(T) vanishes at Tq. 
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We evaluated the equal time correlation of the heat flow by the use of mean-field approximation. This quantity is 
the first term in the Kubo formula and we can obtain a rough estimate for k(T). In the high and low temperature 
limits, our approximation reproduces the result by Harris and Grant 

Similar calculations were also done for the Q2R, a simplified variant of the CCA. The results obtained arc almost 
the same as in the CCA. The normal thermal conductivity was found and it was continuous and smooth at the critical 
point. This proves that the existence of the momentum terms is not relevant to the normal thermal conductivity. 
On the other hand, the dimensionality is important. Energy transport is ballistic in the one-dimensional Q2R. Such 
importance of the dimensionality was reported also in |l3| . 

The similarity of the thermal conductivities in the CCA and the Q2R also implies that the smooth change of n at 
the critical point is rather generic. To investigate to what extent this behavior is generic, however, we must examine 
other dynamical systems with a critical point. It is a future problem. 
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APPENDIX A: 



In this and the next Appendices, we derive the Kubo formula (4.1) for the CCA. We denote a state of the total 
system by to = where w^j = (ciji&ij), and the transformation from the state at time t, u> , to that at time 

t + 1, w t=x , by ft as 

uj t+1 = (Ai) 

Then, the time evolution of any function F(u>) is represented by 

F t+1 (iu) =F*(ft(w)) (A2) 
and F (u>) = F(u>), and the time evolution of a measure p(u>) by 

p<+V)=£*(w,<V))pV) (A3) 



and p (lu) — p{lu), where 



, i 1 if u> = ui' , . , n 

<W)H (A4) 



Now we define the total flux Jtot,x(u) by 

where Jij,a(w) (a = x or y)is the a component of the energy flux at site when the system is in state u>. We 
assume that the initial measure p° equals the local equilibrium measure p\ c defined by 

PleM = in e " AEi ' iM > (^6) 

where E^j(ui) is the local energy around site in state lu. Z\ e denotes the normalization constant 

Zie = J2I[ e ^ iEiAU;) - (A7) 

" hi 







The parameter /3j is the local inverse temperature at the ith column. We consider the temperature variation in the x 
direction only. If the temperature is uniform and all the fts equal a value /3, p\ c becomes the equilibrium measure of 
temperature T = (3" 1 . The average of function F(u>) with respect to the local equilibrium measure is written as 

(F)le = X>(w)Ple(w) (A8) 

Similarly we write the equilibrium average as (F) cq . 

In the following, we calculate the local equilibrium average of the total flux at time t 

t-i 

(«^tot,x)le = ("Aot,x)le + (Jtot.x ~ < Aot,a;)lc (A9) 
t'=0 
t-1 

= (Jtot,*>ie + E £(JfixM - 4*A u ))p» ( A1 °) 

t'=0 w 
t-1 

= (Jtot,x)le + EE Jtot,xH(^H-P°H) (AH) 



t'=0 w 

In the last equality, we have used the identity 

E^ + wm = EE ft+1 H 5 ( w ' = E^^vh- (A12) 

Utilizing the equation of continuity 

^,m(^(<^)) = E^ m (uj) — Jl+l,m, x (uj) + Ji, m ,x(w) ~ Jl,m+l,y(v) + Jl,m,y{^), (A13) 

p 1 ^) is calculated as 

P 1 M = ^E^ fi ( w, ))II e ^ !E! ' m(w ' ) 

uj' l.rn 

x IJexp(-#[i5 ijm (fi(fc/)) + Jj+i.m,^') - ^m,^') + ^,m+i, y (^') - ^m^OD 

/ . in 

= PIc {lo) e n e - A [ j! + i --^')- j! --^')] 

= p le (w) E n(w')) n eW'-P'-M*"^") (A14) 



Inserting the above formula into Eq. ( |All| ), we have 

(4ot,x)le = (Jtot.x)^ + EE4«( W KH | E 5 ( W ' ( w/ ))n e(ft " A " )Jl ' m " (W ' ) - 1 l (A15) 
t'=0 w I u' l,m J 

Now we formally expand the right hand side with respect to VT and obtain in O(VT) 

t-i 

(4ot,x)le - (^tot,x)le + E E E J tot,x(^)Peq(w)5(a;, fi(o/)) ~ 
t-1 

* (Jtot,x)ie + E E J tot;VK q (^) E(# - a-i)^,«,x(«') 

t'=0 i,m 
VT .*i / 

— ("7tot,x)le 7p2 E (^ tot ^^tot,x)oq (A16) 
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where we have used the time invariance of the equilibrium measure, /9 eq (f2(w)) = p cq (ui). 

As we will show in Appendix B, the following equality holds for the first term in the right hand side 

(^tot,x)le — —lyp2\(Jtot,x) )eq ( A 17) 

in O(VT). In addition, we assume that the average energy flux goes to a stationary value in the limit t — > oo 
irrespective of an initial measure. Then the stationary energy flux per site obeys the Fourier heat law 

J st = lim i(4 ot , a )ie = -kVT in O(VT) (A18) 
and the thermal conductivity n is given by 

= E ^v< J to t) x^ot,x)(l - 5<M: (A19) 
t=0 

which is the Kubo formula for the CCA. We remark that the expansion is formal and not justified. The coefficient 
k might be divergent. Currently we have no means to judge the convergence of the coefficient except the numerical 
methods. 



APPENDIX B: 



In this Appendix we prove the formula (A.17). First we note that 
d 



dV/3 



V/3=0 ij !>m 



This is obtained by a straightforward calculation. From now on, we only deal with equilibrium averages and suffix 
"eq" will be omitted. 
Since the total Hamiltonian H is invariant, namely H(Q(uj)) — H(u>), we have 

(Ji,j,xEi, m ) = z- x J2 JijM E i,M e ' mw) 

= Z- iy £j i ^(w)E l , m (u>)e-W a W 

= Z-^^A^i^EiM^ 1 ^))^ "^ (B2) 

where VL^ 1 is the inverse operator of O. Denoting the operation of updating the even sites by O e and that for the odd 
sites by f2°, we can decompose the time evolution operator Q as 

n = n° o n e . (B3) 

Since we have f2 c o f2 c = f2° o f2° = 1 (identity), the inverse operator is given as 

fT 1 = STofT. (B4) 

Let us define the shift operator S by 

(Su)i t j = Wi, 3 -i. (B5) 

This means that the operator 5* shifts the state by one site in the y direction. Because the shift exchanges the roles 
of the even and odd sites, we have 

So!l c = !!%S (B6) 

Sofl° = fl c oS (B7) 

Thus the inverse operator has another representation 
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rr 1 = s^ 1 o n o s 

Inserting this formula into Eq. (B2) and utilizing the shift invariance of the Hamiltonian (i.e., H(S(u>)) 
can write 

(Ji,j, x E ltm ) = Z- 1 ]T JijA^i^El.miS' 1 o n o S(u))e~ pH ^ 

= z- 1 Y JijAn- 1 ° s- x (u))Ei, m +i(n(u))e- f,aM 

From the definition of the flux Ji t j tX (u>) 

j f \ / ^i-i,j( a 'i,j - °m) if (iJ) is even 

iJ ' x[ ] - \ ^i-x,M,3 ~ ^) if (hj) ^ odd, 

where to' = lo = {tOij = {oi.j, o'i.j)}, an d = {^i.j = {<?[ j, <Jj j)}, and the identity SI -1 o = 

we obtain 



in both the cases that site is even or odd. Combining Eqs. (Bll) and (B13), we have 

(Ji,j,x E l,m) = — (Ji t j+i tX E lim+1 ). 



(B8) 
= H(uj)), we 

(B9) 
(BIO) 
(Bll) 

(B12) 

(B13) 
(B14) 



Similarly we have (Ji.j iX ) = ~(Ji.j+i.x) ■ Inserting these equalities into Eq. (Bl), we arrive at 
d 



dV(3 



(Jtot,x)l 



V/3=0 



i,j l.rn 



i,j l,m 

2 ^2 y^J( J hj,x( E lm - E l,m)) 
i,j l,m 

,m,x ~f~ Jl,m,y ^,m+l,y)) 
2" ^ ^ ^ ^ I {Ji,j,x{Jl,m.x ~ ^+l,m,a:)) 
— y ^ ^ ^ (Jj,j,xJl,m,x) 



i,j l,m 
2 ((^tot,x) 2 )- 



(B15) 



This is equivalent to the formula ( |A17| ). 
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FIG. 1. The structure of the system. The two lines at the edges are assigned for heat reservoirs. 

FIG. 2. (a) Distribution of local kinetic energies in Case A. Probability that kinetic energy at a site with horizontal position 
i takes a value is plotted against the value. Calculations were done in a system of size 100 x 100. 
(b) Distribution of kinetic energy in Case B obtained in the same manner as (a) . 

FIG. 3. (a) Scaled temperature profiles in Case A with various system sizes: 30 x 30, 200 x 200, and 300 x 300. 
(b) Scaled temperature profiles for Case B. 



FIG. 4. J t ot,x/L measured in the system of size L x L for various boundary temperatures. The numbers in the figure a-b 
means that Tl = a and Tr — b. 



FIG. 5. The partial Kubo sum k* at various temperatures. 

FIG. 6. Themal conductivity k(T) measured in the direct simulation and that calculated via the Kubo formula in the 
system of size 200 x 200. 



FIG. 7. Thermal conductivity near the critical temperature in the systems with different sizes. 



FIG. 8. Numerically computed k (T) and the mean-field results. 



FIG. 9. Thermal conductivity and 3.5 x k (T). Both arc numerically obtained. 



FIG. 10. Profile of the local energies in the Q2R of various sizes. 



FIG. 11. System size dependence of the total energy flux in the Q2R. 

FIG. 12. Thermal conductivity near the critical point in the Q2R computed via the Kubo formula in the system of size 
100 x 100. 



13 



L-l 



L-2 



C D Q Q Q 



Q Q Q CD 
0-0-0-0 
0-0-0-0 



y direction 



, 



(0000000 0-00000) 

(0000000 (000000) 



o o o o 



L-l 



x direction 



Reservoir 



Fig.l, K. Saito, S. Takesue, S. Miyashita 



14 




2 4 6 8 10 12 

Energy 



Fig. 2a, K. Saito, S. Takesue, S. Miyashita 



15 




Fig. 2b, K. Saito, S. Takesue, S. Miyashita 
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Fig. 3b, K. Saito, S. Takesue, S. Miyashita 
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Fig. 5, K. Saito, S. Takesue, S. Miyashita 
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Fig. 7, K. Saito, S. Takesue, S. Miyashita 
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Fig. 10, K. Saito, S. Takesue, S. Miyashita 
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